/*  sdfdreci.cpp  by K.Tsuru  */
// function ID = 327  DRADIX
/******************************************************
SDouble class 1/a
Provides a reciprocal of a.
[Algorithm]
Get a solusion of an equation a - 1/x = 0 by the
Newton's method with "doubling effective figures"(DEF).
It yileds a recurrence formula
    d=1-a*x(n)
    x(n+1)=x(n)+x(n)*d   2nd
    x(n+1)=x(n)(1+d+d*d) 3rd
    x(n+1)=x(n)(1+d)*(1+d*d) 4th
********************************************************/
#ifndef SN_H
#include "sn.h"
#endif

static const char* const func = "DReciprocal";

#define USE_4TH_CONV_DRec 1

SDouble DReciprocal(const SDouble& x){
//Does not accept BRADIX.
    if( x.Type() != x.REAL ) x.SetError(x.RADIX_ERR, func, 327);
    if( x.Sign(327) == 0 ) x.SetError(x.DIVIDED_BY_ZERO, func, 327);

//Check it can be expressed by a reciprocal of a short number including |x|==1.
    long L, dec_e;
    uint xf = x.Last() - x.First();   // 0.0001 2345 6000(xf = 2)=123456*10^(-9) Ok
    if( (xf <= 2u) &&  x.ConvTolongExp(&L, &dec_e) ){
        SDouble Ds;
        int v = x.Sign(327) > 0 ? 1 : -1;
        Ds = v;
        if( labs(L) != 1 ) Ds = DsDiv(Ds, (ulong)labs(L));      // 1/L, |L| <= SlOpMaxValue()
        if(dec_e) Ds.MultPow10(-dec_e);
        Ds.iterationCount = 0;
        return Ds;
    }
//|x| != 1 in below.
    RealSize C;
    uint max_sz = x.MaxSize();  //not "CurrentMaxSize()"
//delta is an additional value and vx = |x|.
//It takes into irreducible maximum size mode.
    SDouble  y(x.Type(), max_sz), delta(x.Type(), max_sz), vx(Dabs(x));
    int e = vx.NetRdxExp();
//c:counter of iteration, itrmax:its maximum value
    uint c = 0, itrmax = howpow2( (DFIGURES*max_sz)/DOUBLE_FIG+1 ) + 7;
    uint ef = (DOUBLE_FIG*2u)/DFIGURES, fig = C.EffFigures();
//cout << "fig = " << fig << endl;
//Get an approximated value for 1/x.
    vx.MultPowRdx(1-e); //Set exponent at one. 1.0 < vx <= DRADIX
    double dblX = doubleD(vx);
    delta = dblX - vx;
    if( delta.Sign(327) ){
        //If dblX is very close to double value, high precision is necessary.
        ef = max( ef, 2u*(uint)abs(delta.NetRdxExp()) );
    }
    bool fullPrec = false;  //It has evaluated in full perecision or not.
    y = 1.0/dblX;
    y.FixedPoint(y.RdxExp());
    if(ef > fig) ef = fig;
#if USE_4TH_CONV_DRec //==============================
    int pre_exD = 0, now_exD = 0;
//enter into fixed point mode
    do {
        if( (ef = C.SetEffFig(ef)) >= fig) fullPrec = true;
        delta = ONE - vx*y; //evaluate refering CurrentMaxSize()
        now_exD = delta.NetRdxExp();
// In the good converge process, it satisfies now_exD < pre_exD < 0 
        if(now_exD < pre_exD) pre_exD = now_exD;
        else break;     // perhaps comverge

        if(fullPrec && delta.IsMLT(ONE)) {
            C.SetEffFig(0); break;
        }

        y = y*(ONE + delta)*(ONE + delta*delta) ;
        c++;
        ef *= 4u;
        if(ef > fig) ef = fig;
        C.SetEffFig(0);
    } while( (c < itrmax) || !fullPrec );
    //The condition "!delta.isHidden()" corresponds to refering relative error.
    //The condition "c < itrmax" is not necessary but to avoid endless running.
#else // end of USE_4TH_CONV==============================
    // below 2nd convergence
    do {
        if( (ef = C.SetEffFig(ef)) >= fig) fullPrec = true;
        delta = ONE - vx*y; //evaluate refering CurrentMaxSize()
        delta = delta*y;
        y += delta;
        c++;
        ef *= 2u;
        if(ef > fig) ef = fig;
        C.SetEffFig(0);
    } while( ( !delta.IsHidden() && (c < itrmax) ) || !fullPrec );
    //The condition "!delta.isHidden()" corresponds to refering ralative error.
    //The condition "c < itrmax" is not necessary but to avoid endless running.
#endif
    y.PointFree();
    if(c == itrmax) y.SetError(y.FATAL, func, 327);
/*************************************************************
Add a correction.
Let d = 1-x*y'
If d != 0
 y' = (1-d)/x = y*(1-d)  ( y = 1/x )
 y = y'/(1-d) ~= y'*(1+d).
Correction is y'*d
It might not need.
***************************************************************/
#if 1   // Changed from 1 to 0 on Jan 29, 2001.
    delta = ONE - vx*y;
    if(delta.Sign()){   //delta == 0 in almost case.
        ef = delta.Last()-delta.First() +delta.Hidden();
        SDouble ys = y.TakeOutFigures(ef);
        y = y + ys*delta;
    }
#endif
    if(y.Verify()){ // y ~= 1.0
        delta = ONE - y*vx;
        if(!delta.IsMLT(ONE)){
            vx.SetError(vx.VERIFY, func, 327);
        }
 // for debug
 // cout << "in Verify() delta= " << delta << endl;  Wait();
    }
/******************************************************************************
It rounds off if possible.
If this procedure does not exist,a value such as 1/(BRDDIX^2) has an inaccuracy
in the last few figures, though it is a finite decimal. Above correction
cannot remove this inaccuracy.
********************************************************************************/
    if(y.IsPossibleToRound(3)) y.Round(0);
    y.iterationCount = c;
    if(x.Sign() < 0) y.ChangeSign();
    y.MultPowRdx(1-e);      // includes Reform() function.
    return y;
}